##################################################################

# Gautam Nair
# gautam.nair@yale.edu
# Misperceptions of Relative Affluence and Support for International Transfers
# Make Tables 15-20: Subgroup Analysis 
# 15-17 Information vs Control
# 18-20 Information vs Salience

##################################################################

# setwd("")

##################################################################
# loading packages
##################################################################

library(MASS)
library(Hmisc)
library(foreign)
library(sandwich)
library(lmtest)
library(numDeriv)
library(stargazer)
library(ggplot2)
library(plyr)
library(gridExtra)
library(ri)
library(dplyr)
library(scales)

##################################################################

rm(list=ls())

data.working <- readRDS("d_r_cleaned_analysis_dataset.Rda")

vars <- c(
"effort.luck.mix.mostly.luck",
"charity.pretreat.yes",
"identity.cosmopolitan",
"income.50k.atabove",
"ideology.liberal",
"partyid.democrat",
"religion.attendance.often",
"educat.ba.plus",
"work.employed"
)

vars.labels <- c(
"Success is mix of luck and effort or mostly luck",
"Donated to charity in previous year",
"Sees herself as world citizen",
"Household Income>=$50k",
"Liberal",
"Democratic Party ID",
"Often attends religious services",
"Has BA education or higher",
"Employed full time"
)

table.titles <- c("Heterogeneous Treatment Effects on Percentage of Bonus Donated to International Charity (Information vs. Control)",
"Heterogeneous Treatment Effects on Foreign Aid Budget Should be Increased (Information vs. Control)",
"Heterogeneous Treatment Effects on Agricultural Tariffs and Subsidies should be Reduced (Information vs. Control)",
"Heterogeneous Treatment Effects on Percentage of Bonus Donated to International Charity (Information vs. Salience)",
"Heterogeneous Treatment Effects on Foreign Aid Budget Should be Increased (Information vs. Salience)",
"Heterogeneous Treatment Effects on Agricultural Tariffs and Subsidies should be Reduced (Information vs. Salience)"
)

filenames <- c("tf_t_15_subgroup_charity_info_control.txt",
"tf_t_16_subgroup_aid_info_control.txt",
"tf_t_17_subgroup_agro_info_control.txt",
"tf_t_18_subgroup_charity_info_control.txt",
"tf_t_19_subgroup_aid_info_control.txt",
"tf_t_20_subgroup_agro_info_control.txt"
)

first.row.1 <- c("Dummy Variable taking values 1/0",
"CATE (1)",
"CATE (0)",
"Diff. p-val.",
"Control Mean (1)",
"Control Mean (0)",
"N (1)",
"N (0)")

first.row.2 <- c("Dummy Variable taking values 1/0",
"CATE (1)",
"CATE (0)",
"Diff. p-val.",
"Control Mean (1)",
"Control Mean (0)",
"N (1)",
"N (0)")

first.row <- list(first.row.1, first.row.1, first.row.1, first.row.2, first.row.2, first.row.2)

dep.var <- c("q13ab.3", "q8.toolittle", "q11.yes.nodk", "q13ab.3", "q8.toolittle", "q11.yes.nodk")

treatvar <- c("group3")
varslength <- length(vars)
decimals <- 2

temp.tables.list <-vector("list", length(dep.var)) 

z=0

for(y in 1:length(dep.var)){
	z=z+1
	if(z==1|z==2|z==3){
		data.subset <- data.working[data.working$group1==1|data.working$group3==1,]
	}
	if(z==4|z==5|z==6){
		data.subset <- data.working[data.working$group2==1|data.working$group3==1,]
	}
	control.mean <- control.mean.se <- cates <- cate.se <- cate.se.stars <-  pvalue <-  N <- p.val.f <- rep(NA, varslength)
	for (x in c(1, 0)){
		for (i in 1:length(vars)){ 
			index = which(data.subset[,vars[i]]==x)
			tempdata = data.subset[index,]
	
			my.formula <-paste(dep.var[y],'~', paste(treatvar, collapse= ' + '))
			reg.model <- lm(my.formula, data=tempdata)
			cov <- vcovHC(reg.model, type = "HC")
			robustse <- sqrt(diag(cov))
	
			control.mean[i] <- format(round(reg.model$coefficients[1], decimals), nsmall=decimals)
			controlse <- format(round(robustse[1], decimals), nsmall=decimals)
			
			my.formula <-paste(dep.var[y],'~', paste(treatvar, collapse= ' + '))
			ttest <- t.test(tempdata[,dep.var[y]]~tempdata[,treatvar], data=tempdata)
			pvalue[i] <- format(round(ttest$p.value, decimals), nsmall=decimals)
		
			cates[i] <- reg.model$coefficients[2]
			cate.se <- robustse[2]
			cate <- reg.model$coefficients[2]
			cate <- format(round(cate, decimals), nsmall=decimals)
			catese <- format(round(robustse[2], decimals), nsmall=decimals)
	
			cate.se.stars[i] <- paste(cate," ", "(", catese, ")", sep="")
			
			N[i] <- nrow(tempdata)
			
			# F-test of difference
			data.temp <- cbind(data.subset[,treatvar], data.subset[,dep.var[y]], data.subset[,vars[i]])
			Z <- data.temp[,1]
			Y <- data.temp[,2]
			X <- data.temp[,3]
			
			f.pvalue <- anova(lm(Y~X*Z))[3,5]
			p.val.f[i] <-format(round(f.pvalue, decimals), nsmall=decimals)
		}
		if (x==1) {
			cates.1<- cbind (cate.se.stars)
			control.mean.1 <- control.mean
			N.1 <- N
		}
		if (x==0) {
			cates.0<- cbind (cate.se.stars)
			control.mean.0 <- control.mean
			N.0 <-N
		}
	}
	
	temp.table <- cbind(vars.labels, cates.1, cates.0, p.val.f, control.mean.1, control.mean.0, N.1, N.0)
	temp.table <- as.data.frame(temp.table)
	colnames(temp.table) <- first.row[[z]]
	temp.tables.list[[z]] <- temp.table
	stargazer(temp.table, 
	title = table.titles[z],
	out = filenames[z],
	type= "text",
	summary=FALSE, 
	rownames=FALSE)
}	




